Introduction

Short description of the problem

We consider a vacuum vessel embodied as a toroidal volume defined by a 2D arbitrary polygon representing in cross-section in (R,Z) coordinates. We want to efficiently mesh both the volume and the 3D surface of the vessel.

As input, we provide the cross-section polygon and a list dV defining the three resolutions dV = [dR,dZ,dRphi], where each resolution is a distance in meters in the corresponding direction (note the use of dRphi instead of dphi for that purpose). Indeed, we do not want a simple cylindrical mesh, as we want to capture details in the toroidal (Rphi) direction with a similar accuracy for every value of R.

In most uses, we will only need a fraction of the total mesh. Hence, the routine should not compute the total mesh but only the fraction we need. For a given resolution dV (i.e. a given tuple of dR,dZ,dRphi), there should be only one possible mesh. This way, the fraction of interest can be defined either by its physical limits (DR,DZ,DPhi) or by the indices of the mesh elements which are required.

Hence, we need a methods to

Compute the mesh fraction from its limits and return it with its indices
Compute the mesh fraction from its indices and return it with its limits

Not returning the whole mesh makes sense since it may well be a huge array (memory-expensive). The computation should be as fast as possible.

Numerical strategy

Volume

For a given input cross-section polygon VPoly, provided with the required resolutions dR, dZ and dRphi, the total mesh is a torus of rectangular cross-section (R goes from a minimum to a maximum, as Z, and phi goes from -pi to pi). This simplifies the computation of the indices without making the result heavier since only the relevant mesh elements are returned. The limits of the desired submesh can be provided as DR, DZ and DPhi (note that the toroidal limit is angular as opposed to the toroidal resolution).

The minimum and maximim boundaries of R and Z are also provided. They could be deduced from VPOly but they are explicitly provided for two reasons:

In ToFu these limits are stored in the Ves object, so it is already available and there is no need to compute it again
It may happen that a VPoly is slightly modified with no desire to modify the mesh, so it should remain possible to force the boundaries.

The indices are return in numeric form (not boolean), and computed as follows:

The numbering occurs anti-clockwise from -pi to pi for phi
Then, it occurs along Z, from its minimum to its maximum value
Finally it is incremented along R, from its minimum to its maximum value

Hence, we first describe the innermost column (starting from the lowest horizontal ring), and we gradually moved towards the outermost column.

This is necessary since due the desired quasi-constant toroidal accuracy, the number of mesh elements in the toroidal direction is a function of R.

Surface

The strategy for meshing the surface is similar, expect that here there are only two resolutions, dL and dRPhi, where dL is the resolution in the poloidal direction along the segments of VPoly.

Volume meshing

Naive python



In [4]:

    
A = np.random.random((100,))
%timeit len(A)
%timeit A.size









    



The slowest run took 44.64 times longer than the fastest. This could mean that an intermediate result is being cached.
10000000 loops, best of 3: 90.8 ns per loop
The slowest run took 36.53 times longer than the fastest. This could mean that an intermediate result is being cached.
10000000 loops, best of 3: 111 ns per loop

Numpy



In [55]:

    
10%2









    Out[55]:





0



In [201]:

    
import numpy as np

# Preliminary function to get optimal resolution from input resolution
def _Ves_mesh_dlfromL_numpy(LMinMax, dL, DL=None, margin=1.e-9):
    """ Get the actual reolution from the desired resolution and MinMax and limits """
    # Get the number of mesh elements in LMinMax
    N = np.ceil(np.diff(LMinMax)/dL)[0]
    # Derive the real (effective) resolution
    dLr = (np.diff(LMinMax)/N)[0]
    # Get desired limits if any
    DL = LMinMax if DL is None else DL
    # Get the extreme indices of the mesh elements that really need to be created within those limits
    abs0 = np.abs(DL[0]-LMinMax[0])
    if abs0-dLr*np.floor(abs0/dLr)<margin*dLr:
        nL0 = np.round((DL[0]-LMinMax[0])/dLr)
    else:
        nL0 = np.floor((DL[0]-LMinMax[0])/dLr)
    abs1 = np.abs(DL[1]-LMinMax[0])
    if abs1-dLr*np.floor(abs1/dLr)<margin*dLr:
        nL1 = np.round((DL[1]-LMinMax[0])/dLr)-1
    else:
        nL1 = np.floor((DL[1]-LMinMax[0])/dLr)
    # Get the corresponding indices
    indL = np.arange(nL0,nL1+1,1).astype(int)
    # Get the centers of the mesh elements
    L = LMinMax[0] + (0.5 + indL)*dLr
    return L, dLr, indL, N


def _Ves_Vmesh_Phi_numpy(DPhi, dPhir, NRPhi, margin=1.e-9):
    """ Get the min and max indices of the relevant Phi for each R (to be finished) """
    nPhi0 = np.nan*np.ones((dPhir.size,))
    nPhi1 = np.nan*np.ones((dPhir.size,))
    abs0 = np.abs(DPhi[0]+np.pi)
    ind0 = abs0-dPhir*np.floor(abs0/dPhir)<margin*dPhir
    nPhi0[ind0] = np.round((DPhi[0]+np.pi)/dPhir[ind0])
    nPhi0[~ind0] = np.floor((DPhi[0]+np.pi)/dPhir[~ind0])
    abs1 = np.abs(DPhi[1]+np.pi)
    ind1 = abs1-dPhir*np.floor(abs1/dPhir)<margin*dPhir
    nPhi1[ind1] = np.round((DPhi[1]+np.pi)/dPhir[ind1])-1
    nPhi1[~ind1] = np.floor((DPhi[1]+np.pi)/dPhir[~ind1])
    assert np.all(nPhi0>=0) and np.all(nPhi0<=NRPhi)
    assert np.all(nPhi1>=0) and np.all(nPhi1<=NRPhi)
    return nPhi0, nPhi1


def _Ves_Vmesh_Tor_SubFromD_numpy(dR, dZ, dRPhi, RMinMax, ZMinMax, DR=None, DZ=None, DPhi=None, VPoly=None, Out='(X,Y,Z)', margin=1.e-9):
    """ Return the desired submesh indicated by the limits (DR,DZ,DPhi), for the desired resolution (dR,dZ,dRphi) """
    # Get the actual R and Z resolutions and mesh elements
    R0, dRr0, indR0, NR0 = _Ves_mesh_dlfromL_numpy(RMinMax, dR, None, margin=margin)
    R, dRr, indR, NR = _Ves_mesh_dlfromL_numpy(RMinMax, dR, DR, margin=margin)
    Z, dZr, indZ, NZ = _Ves_mesh_dlfromL_numpy(ZMinMax, dZ, DZ, margin=margin)
    Rn, Zn = len(R), len(Z)

    # Get the actual RPhi resolution and Phi mesh elements (! depends on R !)
    NRPhi0 = np.ceil(2.*np.pi*R0/dRPhi)
    NRPhi = np.ceil(2.*np.pi*R/dRPhi)
    dPhir = 2.*np.pi/NRPhi
    # Get the limits if any (and make sure to replace them in the proper quadrants)
    DPhi = [-np.pi,np.pi] if DPhi is None else np.arctan2(np.sin(DPhi),np.cos(DPhi))
    # Later we'll have to remember that the lower limit may be greater than the upper limit due to 2pi modulo
    # Get the two extreme indices of the Phi elements inside the limits for each R
    nPhi0, nPhi1 = _Ves_Vmesh_Phi_numpy(DPhi, dPhir, NRPhi, margin=margin)

    # Compute the centers and indices, using a for loop (no choice) on R
    NRZPhi_cum0 = np.concatenate(([0],NZ*np.cumsum(NRPhi0[:-1])))
    PtsRZP, dV, ind = [], [], []
    for ii in range(0,Rn):
        # Get the phi values
        if DPhi[0]<DPhi[1]:
            indPhi = np.arange(nPhi0[ii],nPhi1[ii]+1)
        else:
            indPhi = np.concatenate((np.arange(nPhi0[ii],NRPhi[ii]),np.arange(0,nPhi1[ii]+1)))
        Phin = len(indPhi)
        phi = -np.pi + (0.5 + indPhi)*dPhir[ii]
        # Get the (R,Z,phi) points
        PtsRZP.append( np.array([R[ii]*np.ones((Phin*Zn,)), np.repeat(Z,Phin), np.tile(phi,Zn)]) )
        # Get the volume
        dV.append( dRr*dZr*dPhir[ii]*R[ii]*np.ones((Phin*Zn,))  )
        # Get the indices
        indii = R0==R[ii]
        ind.append( NRZPhi_cum0[indii] + (indZ[:,np.newaxis]*NRPhi[ii] + indPhi[np.newaxis,:]).flatten() )
        
    PtsRZP = np.concatenate(tuple(PtsRZP),axis=1)
    dV = np.concatenate(tuple(dV))
    ind = np.concatenate(tuple(ind))
    
    # To be commented out once debugged
    assert ind.size==PtsRZP.shape[1]
    assert ind.size==dV.size
    assert np.unique(ind).size==ind.size 

    dRPhir = (dPhir*R)
    if not VPoly is None:
        indin = _Ves_isInside(VPoly, 'Tor', None, PtsRZP[:-1,:], In='(R,Z)')
        PtsRZP, dV, ind = PtsRZP[:,indin], dV[indin], ind[indin]
        Ru = np.unique(PtsRZP[0,:])
        if not Ru==R:
            indii = np.array([R[ii] in Ru for ii in range(0,len(R))], dtype=bool)
            dRPhir = dRPhir[indii]

    if Out.lower()=='(x,y,z)':
        Pts = np.array([PtsRZP[0,:]*np.cos(PtsRZP[2,:]), PtsRZP[0,:]*np.sin(PtsRZP[2,:]), PtsRZP[1,:]])
    else:
        Pts = PtsRZP
        
    return Pts, dV, ind, dRr, dZr, dRPhir

def _Ves_Vmesh_Tor_SubFromInd_numpy(dR, dZ, dRPhi, RMinMax, ZMinMax, ind, Out='(X,Y,Z)', margin=1.e-9):
    """ Return the desired submesh indicated by the (numerical) indices, for the desired resolution (dR,dZ,dRphi) """
    # Get the actual R and Z resolutions and mesh elements
    R, dRr, indR, NR = _Ves_mesh_dlfromL_numpy(RMinMax, dR, None, margin=margin)
    Z, dZr, indZ, NZ = _Ves_mesh_dlfromL_numpy(ZMinMax, dZ, None, margin=margin)
    Rn, Zn = len(R), len(Z)

    # Get the actual RPhi resolution and Phi mesh elements (! depends on R !)
    NRPhi = np.ceil(2.*np.pi*R/dRPhi)
    dPhir = 2.*np.pi/NRPhi
    NRZPhi_cum = np.concatenate(([0],NZ*np.cumsum(NRPhi[:-1])))
    
    # Get R, Z, and Phi indices
    Nind = ind.size
    indR = ind[:,np.newaxis]-NRZPhi_cum[np.newaxis,:]<0
    indRtemp = np.tile(np.arange(0,NR),(Nind,1))
    indRtemp[indR] = np.nan
    indR = np.nanmax(indRtemp, axis=1).astype(int)
    indZ = ((ind - NRZPhi_cum[indR]) // NRPhi[indR]).astype(int)
    indPhi = ind - NRZPhi_cum[indR] - indZ*NRPhi[indR]
    
    # Derive R, Z and phi values
    PtsRZP = np.array([R[indR], Z[indZ], np.nan*np.ones((Nind,))])
    dV = np.nan*np.ones((Nind,))
    Ru = np.unique(PtsRZP[0,:])
    dRPhir = np.nan*np.ones((len(Ru),))
    for ii in range(0,len(Ru)):
        indRii = R==Ru[ii]
        indRu = PtsRZP[0,:]==Ru[ii]
        PtsRZP[2,indRu] = -np.pi + (0.5+indPhi[indRu])*dPhir[R==Ru[ii]]
        dRPhir[ii] = dPhir[indRii]*Ru[ii]
        dV[indRu] = dRr*dZr*dRPhir[ii]
        
    if Out.lower()=='(x,y,z)':
        Pts = np.array([PtsRZP[0,:]*np.cos(PtsRZP[2,:]), PtsRZP[0,:]*np.sin(PtsRZP[2,:]), PtsRZP[1,:]])
    else:
        Pts = PtsRZP

    return Pts, dV, dRr, dZr, dRPhir

Check and debug



In [202]:

    
# Check and debug
[dR, dZ, dRPhi], RMinMax, ZMinMax = [0.03]*3, [2.,3.], [-1.,1.]
DR, DZ, DPhi = [2.15,2.6], [-0.21,0.51], [7.*np.pi/4.,5.*np.pi/4.]

Pts0, dV0, ind0, dRr0, dZr0, dRPhir0 = _Ves_Vmesh_Tor_SubFromD_numpy(dR,dZ,dRPhi, RMinMax,ZMinMax, Out='(X,Y,Z)')
PtsRZP0 = np.array([np.hypot(Pts0[0,:],Pts0[1,:]), Pts0[2,:],np.arctan2(Pts0[1,:],Pts0[0,:])])
Pts0b, dV0b, dRr0b, dZr0b, dRPhir0b = _Ves_Vmesh_Tor_SubFromInd_numpy(dR, dZ, dRPhi, RMinMax, ZMinMax, ind0, Out='(X,Y,Z)')
print "Check 0", [np.allclose(Pts0,Pts0b,atol=1.e-14,rtol=0.,equal_nan=True) for (A,B) in [(Pts0,Pts0b),(dV0,dV0b),(dRr0,dRr0b),(dZr0,dZr0b),(dRPhir0,dRPhir0b)]]

Pts, dV, ind, dRr, dZr, dRPhir = _Ves_Vmesh_Tor_SubFromD_numpy(dR,dZ,dRPhi, RMinMax,ZMinMax, DR=DR,DZ=DZ,DPhi=DPhi, Out='(X,Y,Z)')
PtsRZP = np.array([np.hypot(Pts[0,:],Pts[1,:]), Pts[2,:],np.arctan2(Pts[1,:],Pts[0,:])])
Ptsb, dVb, dRrb, dZrb, dRPhirb = _Ves_Vmesh_Tor_SubFromInd_numpy(dR,dZ,dRPhi, RMinMax,ZMinMax, ind, Out='(X,Y,Z)')
print "Check 1", [np.allclose(A,B,atol=1.e-14,rtol=0.,equal_nan=True) for (A,B) in [(Pts,Ptsb),(dV,dVb),(dRr,dRrb),(dZr,dZrb),(dRPhir,dRPhirb)]]









    



Check 0 [True, True, True, True, True]
Check 1 [True, True, True, True, True]



In [203]:

    
# Plot for check and debug
import matplotlib.pyplot as plt

plt.figure(figsize=(16,10))
plt.subplot(1,2,1, aspect="equal", adjustable="datalim")
plt.plot(PtsRZP0[0,:],PtsRZP0[1,:],'.b', PtsRZP[0,:],PtsRZP[1,:],'.r', markersize=2)
plt.plot([DR[0],DR[1],DR[1],DR[0],DR[0]], [DZ[0],DZ[0],DZ[1],DZ[1],DZ[0]],'--k')

plt.subplot(1,2,2, aspect="equal", adjustable="datalim")
plt.plot(Pts0[0,:],Pts0[1,:],'.b', Pts[0,:],Pts[1,:],'.r', markersize=2)
thet = np.linspace(DPhi[0],DPhi[1],100) if DPhi[0]<DPhi[1] else np.linspace(DPhi[0],DPhi[1]+2.*np.pi,100)
X = np.concatenate((DR[0]*np.cos(thet),DR[1]*np.cos(thet[::-1]),[DR[0]*np.cos(thet[0])]))
Y = np.concatenate((DR[0]*np.sin(thet),DR[1]*np.sin(thet[::-1]),[DR[0]*np.sin(thet[0])]))
plt.plot(X,Y,'--k')

plt.show()



In [122]:

    
# Check and debug
[dR, dZ, dRPhi], RMinMax, ZMinMax = [0.02]*3, [2.,3.], [-1.,1.]
DR, DZ, DPhi = [2.15,2.6], [-0.21,0.51], [3.*np.pi/4.,-np.pi/4.]
Pts0, dV0, ind0, dRr0, dZr0, dRPhir0 = _Ves_Vmesh_Tor_SubFromD(dR,dZ,dRPhi, RMinMax,ZMinMax, Out='(X,Y,Z)')

%timeit out = _Ves_Vmesh_Tor_SubFromD(dR,dZ,dRPhi, RMinMax,ZMinMax, Out='(X,Y,Z)')
%timeit out = _Ves_Vmesh_Tor_SubFromInd(dR, dZ, dRPhi, RMinMax, ZMinMax, ind0, Out='(X,Y,Z)')









    



1 loop, best of 3: 606 ms per loop
1 loop, best of 3: 5.07 s per loop

Cython



In [206]:

    
%load_ext cython



In [219]:

    
%%cython -a
cimport cython
import numpy as np
cimport numpy as cnp
from cpython cimport bool
from libc.math cimport ceil as Cceil, abs as Cabs, floor as Cfloor, round as Cround

#@cython.cdivision(True)
#@cython.wraparound(False)
#@cython.boundscheck(False)

# Preliminary function to get optimal resolution from input resolution
@cython.cdivision(True)
@cython.wraparound(False)
@cython.boundscheck(False)
def _Ves_mesh_dlfromL_cython(double[::1] LMinMax, double dL, DL=None, double margin=1.e-9):
    """ Get the actual reolution from the desired resolution and MinMax and limits """
    # Get the number of mesh elements in LMinMax
    cdef double N = Cceil((LMinMax[1]-LMinMax[0])/dL)
    # Derive the real (effective) resolution
    cdef double dLr = (LMinMax[1]-LMinMax[0])/N
    # Get desired limits if any
    cdef double[::1] DLc, L
    cdef int[::1] indL
    cdef double nL0, nL1, abs0, abs1
    if DL is None:
        DLc = LMinMax
    else:
        DLc = DL
    # Get the extreme indices of the mesh elements that really need to be created within those limits
    abs0 = Cabs(DLc[0]-LMinMax[0])
    if abs0-dLr*Cfloor(abs0/dLr)<margin*dLr:
        nL0 = Cround((DLc[0]-LMinMax[0])/dLr)
    else:
        nL0 = Cfloor((DLc[0]-LMinMax[0])/dLr)
    abs1 = Cabs(DLc[1]-LMinMax[0])
    if abs1-dLr*Cfloor(abs1/dLr)<margin*dLr:
        nL1 = Cround((DLc[1]-LMinMax[0])/dLr)-1
    else:
        nL1 = Cfloor((DLc[1]-LMinMax[0])/dLr)
    # Get the corresponding indices
    for ii in range(nL0,nL1+1):
    indL = np.arange(nL0,nL1+1,1).astype(int)
    # Get the centers of the mesh elements
    L = LMinMax[0] + (0.5 + indL)*dLr
    return L, dLr, indL, N


def _Ves_Vmesh_Phi_numpy(DPhi, dPhir, NRPhi, margin=1.e-9):
    """ Get the min and max indices of the relevant Phi for each R (to be finished) """
    nPhi0 = np.nan*np.ones((dPhir.size,))
    nPhi1 = np.nan*np.ones((dPhir.size,))
    abs0 = np.abs(DPhi[0]+np.pi)
    ind0 = abs0-dPhir*np.floor(abs0/dPhir)<margin*dPhir
    nPhi0[ind0] = np.round((DPhi[0]+np.pi)/dPhir[ind0])
    nPhi0[~ind0] = np.floor((DPhi[0]+np.pi)/dPhir[~ind0])
    abs1 = np.abs(DPhi[1]+np.pi)
    ind1 = abs1-dPhir*np.floor(abs1/dPhir)<margin*dPhir
    nPhi1[ind1] = np.round((DPhi[1]+np.pi)/dPhir[ind1])-1
    nPhi1[~ind1] = np.floor((DPhi[1]+np.pi)/dPhir[~ind1])
    assert np.all(nPhi0>=0) and np.all(nPhi0<=NRPhi)
    assert np.all(nPhi1>=0) and np.all(nPhi1<=NRPhi)
    return nPhi0, nPhi1









    Out[219]:









    
    Cython: _cython_magic_1e9ae0b4a9168851b26e8dfa94f3a771.pyx
    
    


Generated by Cython 0.25.2

    Yellow lines hint at Python interaction.

    Click on a line that starts with a "+" to see the C code that Cython generated for it.

+01: cimport cython
  __pyx_t_1 = PyDict_New(); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_test, __pyx_t_1) < 0) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+02: import numpy as np
  __pyx_t_1 = __Pyx_Import(__pyx_n_s_numpy, 0, -1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 2, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_np, __pyx_t_1) < 0) __PYX_ERR(0, 2, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
 03: cimport numpy as cnp
 04: from cpython cimport bool
 05: from libc.math cimport ceil as Cceil, abs as Cabs, floor as Cfloor, round as Cround
 06: 
 07: #@cython.cdivision(True)
 08: #@cython.wraparound(False)
 09: #@cython.boundscheck(False)
 10: 
 11: # Preliminary function to get optimal resolution from input resolution
 12: @cython.cdivision(True)
 13: @cython.wraparound(False)
 14: @cython.boundscheck(False)
+15: def _Ves_mesh_dlfromL_cython(double[::1] LMinMax, double dL, DL=None, double margin=1.e-9):
/* Python wrapper */
static PyObject *__pyx_pw_46_cython_magic_1e9ae0b4a9168851b26e8dfa94f3a771_1_Ves_mesh_dlfromL_cython(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
static char __pyx_doc_46_cython_magic_1e9ae0b4a9168851b26e8dfa94f3a771__Ves_mesh_dlfromL_cython[] = " Get the actual reolution from the desired resolution and MinMax and limits ";
static PyMethodDef __pyx_mdef_46_cython_magic_1e9ae0b4a9168851b26e8dfa94f3a771_1_Ves_mesh_dlfromL_cython = {"_Ves_mesh_dlfromL_cython", (PyCFunction)__pyx_pw_46_cython_magic_1e9ae0b4a9168851b26e8dfa94f3a771_1_Ves_mesh_dlfromL_cython, METH_VARARGS|METH_KEYWORDS, __pyx_doc_46_cython_magic_1e9ae0b4a9168851b26e8dfa94f3a771__Ves_mesh_dlfromL_cython};
static PyObject *__pyx_pw_46_cython_magic_1e9ae0b4a9168851b26e8dfa94f3a771_1_Ves_mesh_dlfromL_cython(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
  __Pyx_memviewslice __pyx_v_LMinMax = { 0, 0, { 0 }, { 0 }, { 0 } };
  double __pyx_v_dL;
  PyObject *__pyx_v_DL = 0;
  double __pyx_v_margin;
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("_Ves_mesh_dlfromL_cython (wrapper)", 0);
  {
    static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_LMinMax,&__pyx_n_s_dL,&__pyx_n_s_DL,&__pyx_n_s_margin,0};
    PyObject* values[4] = {0,0,0,0};
    values[2] = ((PyObject *)Py_None);
    if (unlikely(__pyx_kwds)) {
      Py_ssize_t kw_args;
      const Py_ssize_t pos_args = PyTuple_GET_SIZE(__pyx_args);
      switch (pos_args) {
        case  4: values[3] = PyTuple_GET_ITEM(__pyx_args, 3);
        case  3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
        case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
        case  1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
        case  0: break;
        default: goto __pyx_L5_argtuple_error;
      }
      kw_args = PyDict_Size(__pyx_kwds);
      switch (pos_args) {
        case  0:
        if (likely((values[0] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_LMinMax)) != 0)) kw_args--;
        else goto __pyx_L5_argtuple_error;
        case  1:
        if (likely((values[1] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_dL)) != 0)) kw_args--;
        else {
          __Pyx_RaiseArgtupleInvalid("_Ves_mesh_dlfromL_cython", 0, 2, 4, 1); __PYX_ERR(0, 15, __pyx_L3_error)
        }
        case  2:
        if (kw_args > 0) {
          PyObject* value = PyDict_GetItem(__pyx_kwds, __pyx_n_s_DL);
          if (value) { values[2] = value; kw_args--; }
        }
        case  3:
        if (kw_args > 0) {
          PyObject* value = PyDict_GetItem(__pyx_kwds, __pyx_n_s_margin);
          if (value) { values[3] = value; kw_args--; }
        }
      }
      if (unlikely(kw_args > 0)) {
        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "_Ves_mesh_dlfromL_cython") < 0)) __PYX_ERR(0, 15, __pyx_L3_error)
      }
    } else {
      switch (PyTuple_GET_SIZE(__pyx_args)) {
        case  4: values[3] = PyTuple_GET_ITEM(__pyx_args, 3);
        case  3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
        case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
        values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
        break;
        default: goto __pyx_L5_argtuple_error;
      }
    }
    __pyx_v_LMinMax = __Pyx_PyObject_to_MemoryviewSlice_dc_double(values[0]); if (unlikely(!__pyx_v_LMinMax.memview)) __PYX_ERR(0, 15, __pyx_L3_error)
    __pyx_v_dL = __pyx_PyFloat_AsDouble(values[1]); if (unlikely((__pyx_v_dL == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 15, __pyx_L3_error)
    __pyx_v_DL = values[2];
    if (values[3]) {
      __pyx_v_margin = __pyx_PyFloat_AsDouble(values[3]); if (unlikely((__pyx_v_margin == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 15, __pyx_L3_error)
    } else {
      __pyx_v_margin = ((double)1.e-9);
    }
  }
  goto __pyx_L4_argument_unpacking_done;
  __pyx_L5_argtuple_error:;
  __Pyx_RaiseArgtupleInvalid("_Ves_mesh_dlfromL_cython", 0, 2, 4, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 15, __pyx_L3_error)
  __pyx_L3_error:;
  __Pyx_AddTraceback("_cython_magic_1e9ae0b4a9168851b26e8dfa94f3a771._Ves_mesh_dlfromL_cython", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __Pyx_RefNannyFinishContext();
  return NULL;
  __pyx_L4_argument_unpacking_done:;
  __pyx_r = __pyx_pf_46_cython_magic_1e9ae0b4a9168851b26e8dfa94f3a771__Ves_mesh_dlfromL_cython(__pyx_self, __pyx_v_LMinMax, __pyx_v_dL, __pyx_v_DL, __pyx_v_margin);

  /* function exit code */
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_pf_46_cython_magic_1e9ae0b4a9168851b26e8dfa94f3a771__Ves_mesh_dlfromL_cython(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_LMinMax, double __pyx_v_dL, PyObject *__pyx_v_DL, double __pyx_v_margin) {
  double __pyx_v_N;
  double __pyx_v_dLr;
  __Pyx_memviewslice __pyx_v_DLc = { 0, 0, { 0 }, { 0 }, { 0 } };
  double __pyx_v_nL0;
  double __pyx_v_nL1;
  double __pyx_v_abs0;
  double __pyx_v_abs1;
  PyObject *__pyx_v_indL = NULL;
  PyObject *__pyx_v_L = NULL;
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("_Ves_mesh_dlfromL_cython", 0);
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __PYX_XDEC_MEMVIEW(&__pyx_t_7, 1);
  __Pyx_XDECREF(__pyx_t_20);
  __Pyx_XDECREF(__pyx_t_21);
  __Pyx_XDECREF(__pyx_t_22);
  __Pyx_XDECREF(__pyx_t_23);
  __Pyx_XDECREF(__pyx_t_24);
  __Pyx_XDECREF(__pyx_t_25);
  __Pyx_XDECREF(__pyx_t_27);
  __Pyx_AddTraceback("_cython_magic_1e9ae0b4a9168851b26e8dfa94f3a771._Ves_mesh_dlfromL_cython", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  __pyx_L0:;
  __PYX_XDEC_MEMVIEW(&__pyx_v_DLc, 1);
  __Pyx_XDECREF(__pyx_v_indL);
  __Pyx_XDECREF(__pyx_v_L);
  __PYX_XDEC_MEMVIEW(&__pyx_v_LMinMax, 1);
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
/* … */
  __pyx_tuple__23 = PyTuple_Pack(13, __pyx_n_s_LMinMax, __pyx_n_s_dL, __pyx_n_s_DL, __pyx_n_s_margin, __pyx_n_s_N, __pyx_n_s_dLr, __pyx_n_s_DLc, __pyx_n_s_nL0, __pyx_n_s_nL1, __pyx_n_s_abs0, __pyx_n_s_abs1, __pyx_n_s_indL, __pyx_n_s_L); if (unlikely(!__pyx_tuple__23)) __PYX_ERR(0, 15, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_tuple__23);
  __Pyx_GIVEREF(__pyx_tuple__23);
/* … */
  __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_46_cython_magic_1e9ae0b4a9168851b26e8dfa94f3a771_1_Ves_mesh_dlfromL_cython, NULL, __pyx_n_s_cython_magic_1e9ae0b4a9168851b2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 15, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_Ves_mesh_dlfromL_cython, __pyx_t_1) < 0) __PYX_ERR(0, 15, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_codeobj__24 = (PyObject*)__Pyx_PyCode_New(4, 0, 13, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__23, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_Home_DV226270_cache_ipython_cyt, __pyx_n_s_Ves_mesh_dlfromL_cython, 15, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__24)) __PYX_ERR(0, 15, __pyx_L1_error)
 16:     """ Get the actual reolution from the desired resolution and MinMax and limits """
 17:     # Get the number of mesh elements in LMinMax
+18:     cdef double N = Cceil((LMinMax[1]-LMinMax[0])/dL)
  __pyx_t_1 = 1;
  __pyx_t_2 = 0;
  __pyx_v_N = ceil((((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_LMinMax.data) + __pyx_t_1)) ))) - (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_LMinMax.data) + __pyx_t_2)) )))) / __pyx_v_dL));
 19:     # Derive the real (effective) resolution
+20:     cdef double dLr = (LMinMax[1]-LMinMax[0])/N
  __pyx_t_3 = 1;
  __pyx_t_4 = 0;
  __pyx_v_dLr = (((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_LMinMax.data) + __pyx_t_3)) ))) - (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_LMinMax.data) + __pyx_t_4)) )))) / __pyx_v_N);
 21:     # Get desired limits if any
 22:     cdef double[::1] DLc
 23:     cdef double nL0, nL1, abs0, abs1
+24:     if DL is None:
  __pyx_t_5 = (__pyx_v_DL == Py_None);
  __pyx_t_6 = (__pyx_t_5 != 0);
  if (__pyx_t_6) {
/* … */
    goto __pyx_L3;
  }
+25:         DLc = LMinMax
    __PYX_INC_MEMVIEW(&__pyx_v_LMinMax, 0);
    __pyx_v_DLc = __pyx_v_LMinMax;
 26:     else:
+27:         DLc = DL
  /*else*/ {
    __pyx_t_7 = __Pyx_PyObject_to_MemoryviewSlice_dc_double(__pyx_v_DL);
    if (unlikely(!__pyx_t_7.memview)) __PYX_ERR(0, 27, __pyx_L1_error)
    __pyx_v_DLc = __pyx_t_7;
    __pyx_t_7.memview = NULL;
    __pyx_t_7.data = NULL;
  }
  __pyx_L3:;
 28:     # Get the extreme indices of the mesh elements that really need to be created within those limits
+29:     abs0 = Cabs(DLc[0]-LMinMax[0])
  __pyx_t_8 = 0;
  __pyx_t_9 = 0;
  __pyx_v_abs0 = fabs(((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_DLc.data) + __pyx_t_8)) ))) - (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_LMinMax.data) + __pyx_t_9)) )))));
+30:     if abs0-dLr*Cfloor(abs0/dLr)<margin*dLr:
  __pyx_t_6 = (((__pyx_v_abs0 - (__pyx_v_dLr * floor((__pyx_v_abs0 / __pyx_v_dLr)))) < (__pyx_v_margin * __pyx_v_dLr)) != 0);
  if (__pyx_t_6) {
/* … */
    goto __pyx_L4;
  }
+31:         nL0 = Cround((DLc[0]-LMinMax[0])/dLr)
    __pyx_t_10 = 0;
    __pyx_t_11 = 0;
    __pyx_v_nL0 = round((((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_DLc.data) + __pyx_t_10)) ))) - (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_LMinMax.data) + __pyx_t_11)) )))) / __pyx_v_dLr));
 32:     else:
+33:         nL0 = Cfloor((DLc[0]-LMinMax[0])/dLr)
  /*else*/ {
    __pyx_t_12 = 0;
    __pyx_t_13 = 0;
    __pyx_v_nL0 = floor((((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_DLc.data) + __pyx_t_12)) ))) - (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_LMinMax.data) + __pyx_t_13)) )))) / __pyx_v_dLr));
  }
  __pyx_L4:;
+34:     abs1 = Cabs(DLc[1]-LMinMax[0])
  __pyx_t_14 = 1;
  __pyx_t_15 = 0;
  __pyx_v_abs1 = fabs(((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_DLc.data) + __pyx_t_14)) ))) - (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_LMinMax.data) + __pyx_t_15)) )))));
+35:     if abs1-dLr*Cfloor(abs1/dLr)<margin*dLr:
  __pyx_t_6 = (((__pyx_v_abs1 - (__pyx_v_dLr * floor((__pyx_v_abs1 / __pyx_v_dLr)))) < (__pyx_v_margin * __pyx_v_dLr)) != 0);
  if (__pyx_t_6) {
/* … */
    goto __pyx_L5;
  }
+36:         nL1 = Cround((DLc[1]-LMinMax[0])/dLr)-1
    __pyx_t_16 = 1;
    __pyx_t_17 = 0;
    __pyx_v_nL1 = (round((((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_DLc.data) + __pyx_t_16)) ))) - (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_LMinMax.data) + __pyx_t_17)) )))) / __pyx_v_dLr)) - 1.0);
 37:     else:
+38:         nL1 = Cfloor((DLc[1]-LMinMax[0])/dLr)
  /*else*/ {
    __pyx_t_18 = 1;
    __pyx_t_19 = 0;
    __pyx_v_nL1 = floor((((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_DLc.data) + __pyx_t_18)) ))) - (*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_LMinMax.data) + __pyx_t_19)) )))) / __pyx_v_dLr));
  }
  __pyx_L5:;
 39:     # Get the corresponding indices
+40:     indL = np.arange(nL0,nL1+1,1).astype(int)
  __pyx_t_22 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_22)) __PYX_ERR(0, 40, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_22);
  __pyx_t_23 = __Pyx_PyObject_GetAttrStr(__pyx_t_22, __pyx_n_s_arange); if (unlikely(!__pyx_t_23)) __PYX_ERR(0, 40, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_23);
  __Pyx_DECREF(__pyx_t_22); __pyx_t_22 = 0;
  __pyx_t_22 = PyFloat_FromDouble(__pyx_v_nL0); if (unlikely(!__pyx_t_22)) __PYX_ERR(0, 40, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_22);
  __pyx_t_24 = PyFloat_FromDouble((__pyx_v_nL1 + 1.0)); if (unlikely(!__pyx_t_24)) __PYX_ERR(0, 40, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_24);
  __pyx_t_25 = NULL;
  __pyx_t_26 = 0;
  if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_23))) {
    __pyx_t_25 = PyMethod_GET_SELF(__pyx_t_23);
    if (likely(__pyx_t_25)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_23);
      __Pyx_INCREF(__pyx_t_25);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_23, function);
      __pyx_t_26 = 1;
    }
  }
  #if CYTHON_FAST_PYCALL
  if (PyFunction_Check(__pyx_t_23)) {
    PyObject *__pyx_temp[4] = {__pyx_t_25, __pyx_t_22, __pyx_t_24, __pyx_int_1};
    __pyx_t_21 = __Pyx_PyFunction_FastCall(__pyx_t_23, __pyx_temp+1-__pyx_t_26, 3+__pyx_t_26); if (unlikely(!__pyx_t_21)) __PYX_ERR(0, 40, __pyx_L1_error)
    __Pyx_XDECREF(__pyx_t_25); __pyx_t_25 = 0;
    __Pyx_GOTREF(__pyx_t_21);
    __Pyx_DECREF(__pyx_t_22); __pyx_t_22 = 0;
    __Pyx_DECREF(__pyx_t_24); __pyx_t_24 = 0;
  } else
  #endif
  #if CYTHON_FAST_PYCCALL
  if (__Pyx_PyFastCFunction_Check(__pyx_t_23)) {
    PyObject *__pyx_temp[4] = {__pyx_t_25, __pyx_t_22, __pyx_t_24, __pyx_int_1};
    __pyx_t_21 = __Pyx_PyCFunction_FastCall(__pyx_t_23, __pyx_temp+1-__pyx_t_26, 3+__pyx_t_26); if (unlikely(!__pyx_t_21)) __PYX_ERR(0, 40, __pyx_L1_error)
    __Pyx_XDECREF(__pyx_t_25); __pyx_t_25 = 0;
    __Pyx_GOTREF(__pyx_t_21);
    __Pyx_DECREF(__pyx_t_22); __pyx_t_22 = 0;
    __Pyx_DECREF(__pyx_t_24); __pyx_t_24 = 0;
  } else
  #endif
  {
    __pyx_t_27 = PyTuple_New(3+__pyx_t_26); if (unlikely(!__pyx_t_27)) __PYX_ERR(0, 40, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_27);
    if (__pyx_t_25) {
      __Pyx_GIVEREF(__pyx_t_25); PyTuple_SET_ITEM(__pyx_t_27, 0, __pyx_t_25); __pyx_t_25 = NULL;
    }
    __Pyx_GIVEREF(__pyx_t_22);
    PyTuple_SET_ITEM(__pyx_t_27, 0+__pyx_t_26, __pyx_t_22);
    __Pyx_GIVEREF(__pyx_t_24);
    PyTuple_SET_ITEM(__pyx_t_27, 1+__pyx_t_26, __pyx_t_24);
    __Pyx_INCREF(__pyx_int_1);
    __Pyx_GIVEREF(__pyx_int_1);
    PyTuple_SET_ITEM(__pyx_t_27, 2+__pyx_t_26, __pyx_int_1);
    __pyx_t_22 = 0;
    __pyx_t_24 = 0;
    __pyx_t_21 = __Pyx_PyObject_Call(__pyx_t_23, __pyx_t_27, NULL); if (unlikely(!__pyx_t_21)) __PYX_ERR(0, 40, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_21);
    __Pyx_DECREF(__pyx_t_27); __pyx_t_27 = 0;
  }
  __Pyx_DECREF(__pyx_t_23); __pyx_t_23 = 0;
  __pyx_t_23 = __Pyx_PyObject_GetAttrStr(__pyx_t_21, __pyx_n_s_astype); if (unlikely(!__pyx_t_23)) __PYX_ERR(0, 40, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_23);
  __Pyx_DECREF(__pyx_t_21); __pyx_t_21 = 0;
  __pyx_t_21 = NULL;
  if (CYTHON_UNPACK_METHODS && likely(PyMethod_Check(__pyx_t_23))) {
    __pyx_t_21 = PyMethod_GET_SELF(__pyx_t_23);
    if (likely(__pyx_t_21)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_23);
      __Pyx_INCREF(__pyx_t_21);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_23, function);
    }
  }
  if (!__pyx_t_21) {
    __pyx_t_20 = __Pyx_PyObject_CallOneArg(__pyx_t_23, ((PyObject *)(&PyInt_Type))); if (unlikely(!__pyx_t_20)) __PYX_ERR(0, 40, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_20);
  } else {
    #if CYTHON_FAST_PYCALL
    if (PyFunction_Check(__pyx_t_23)) {
      PyObject *__pyx_temp[2] = {__pyx_t_21, ((PyObject *)(&PyInt_Type))};
      __pyx_t_20 = __Pyx_PyFunction_FastCall(__pyx_t_23, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_20)) __PYX_ERR(0, 40, __pyx_L1_error)
      __Pyx_XDECREF(__pyx_t_21); __pyx_t_21 = 0;
      __Pyx_GOTREF(__pyx_t_20);
    } else
    #endif
    #if CYTHON_FAST_PYCCALL
    if (__Pyx_PyFastCFunction_Check(__pyx_t_23)) {
      PyObject *__pyx_temp[2] = {__pyx_t_21, ((PyObject *)(&PyInt_Type))};
      __pyx_t_20 = __Pyx_PyCFunction_FastCall(__pyx_t_23, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_20)) __PYX_ERR(0, 40, __pyx_L1_error)
      __Pyx_XDECREF(__pyx_t_21); __pyx_t_21 = 0;
      __Pyx_GOTREF(__pyx_t_20);
    } else
    #endif
    {
      __pyx_t_27 = PyTuple_New(1+1); if (unlikely(!__pyx_t_27)) __PYX_ERR(0, 40, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_27);
      __Pyx_GIVEREF(__pyx_t_21); PyTuple_SET_ITEM(__pyx_t_27, 0, __pyx_t_21); __pyx_t_21 = NULL;
      __Pyx_INCREF(((PyObject *)(&PyInt_Type)));
      __Pyx_GIVEREF(((PyObject *)(&PyInt_Type)));
      PyTuple_SET_ITEM(__pyx_t_27, 0+1, ((PyObject *)(&PyInt_Type)));
      __pyx_t_20 = __Pyx_PyObject_Call(__pyx_t_23, __pyx_t_27, NULL); if (unlikely(!__pyx_t_20)) __PYX_ERR(0, 40, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_20);
      __Pyx_DECREF(__pyx_t_27); __pyx_t_27 = 0;
    }
  }
  __Pyx_DECREF(__pyx_t_23); __pyx_t_23 = 0;
  __pyx_v_indL = __pyx_t_20;
  __pyx_t_20 = 0;
 41:     # Get the centers of the mesh elements
+42:     L = LMinMax[0] + (0.5 + indL)*dLr
  __pyx_t_28 = 0;
  __pyx_t_20 = PyFloat_FromDouble((*((double *) ( /* dim=0 */ ((char *) (((double *) __pyx_v_LMinMax.data) + __pyx_t_28)) )))); if (unlikely(!__pyx_t_20)) __PYX_ERR(0, 42, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_20);
  __pyx_t_23 = __Pyx_PyFloat_AddCObj(__pyx_float_0_5, __pyx_v_indL, 0.5, 0); if (unlikely(!__pyx_t_23)) __PYX_ERR(0, 42, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_23);
  __pyx_t_27 = PyFloat_FromDouble(__pyx_v_dLr); if (unlikely(!__pyx_t_27)) __PYX_ERR(0, 42, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_27);
  __pyx_t_21 = PyNumber_Multiply(__pyx_t_23, __pyx_t_27); if (unlikely(!__pyx_t_21)) __PYX_ERR(0, 42, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_21);
  __Pyx_DECREF(__pyx_t_23); __pyx_t_23 = 0;
  __Pyx_DECREF(__pyx_t_27); __pyx_t_27 = 0;
  __pyx_t_27 = PyNumber_Add(__pyx_t_20, __pyx_t_21); if (unlikely(!__pyx_t_27)) __PYX_ERR(0, 42, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_27);
  __Pyx_DECREF(__pyx_t_20); __pyx_t_20 = 0;
  __Pyx_DECREF(__pyx_t_21); __pyx_t_21 = 0;
  __pyx_v_L = __pyx_t_27;
  __pyx_t_27 = 0;
+43:     return L, dLr, indL, N
  __Pyx_XDECREF(__pyx_r);
  __pyx_t_27 = PyFloat_FromDouble(__pyx_v_dLr); if (unlikely(!__pyx_t_27)) __PYX_ERR(0, 43, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_27);
  __pyx_t_21 = PyFloat_FromDouble(__pyx_v_N); if (unlikely(!__pyx_t_21)) __PYX_ERR(0, 43, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_21);
  __pyx_t_20 = PyTuple_New(4); if (unlikely(!__pyx_t_20)) __PYX_ERR(0, 43, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_20);
  __Pyx_INCREF(__pyx_v_L);
  __Pyx_GIVEREF(__pyx_v_L);
  PyTuple_SET_ITEM(__pyx_t_20, 0, __pyx_v_L);
  __Pyx_GIVEREF(__pyx_t_27);
  PyTuple_SET_ITEM(__pyx_t_20, 1, __pyx_t_27);
  __Pyx_INCREF(__pyx_v_indL);
  __Pyx_GIVEREF(__pyx_v_indL);
  PyTuple_SET_ITEM(__pyx_t_20, 2, __pyx_v_indL);
  __Pyx_GIVEREF(__pyx_t_21);
  PyTuple_SET_ITEM(__pyx_t_20, 3, __pyx_t_21);
  __pyx_t_27 = 0;
  __pyx_t_21 = 0;
  __pyx_r = __pyx_t_20;
  __pyx_t_20 = 0;
  goto __pyx_L0;
 44: 
 45: 
+46: def _Ves_Vmesh_Phi_numpy(DPhi, dPhir, NRPhi, margin=1.e-9):
/* Python wrapper */
static PyObject *__pyx_pw_46_cython_magic_1e9ae0b4a9168851b26e8dfa94f3a771_3_Ves_Vmesh_Phi_numpy(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
static char __pyx_doc_46_cython_magic_1e9ae0b4a9168851b26e8dfa94f3a771_2_Ves_Vmesh_Phi_numpy[] = " Get the min and max indices of the relevant Phi for each R (to be finished) ";
static PyMethodDef __pyx_mdef_46_cython_magic_1e9ae0b4a9168851b26e8dfa94f3a771_3_Ves_Vmesh_Phi_numpy = {"_Ves_Vmesh_Phi_numpy", (PyCFunction)__pyx_pw_46_cython_magic_1e9ae0b4a9168851b26e8dfa94f3a771_3_Ves_Vmesh_Phi_numpy, METH_VARARGS|METH_KEYWORDS, __pyx_doc_46_cython_magic_1e9ae0b4a9168851b26e8dfa94f3a771_2_Ves_Vmesh_Phi_numpy};
static PyObject *__pyx_pw_46_cython_magic_1e9ae0b4a9168851b26e8dfa94f3a771_3_Ves_Vmesh_Phi_numpy(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
  PyObject *__pyx_v_DPhi = 0;
  PyObject *__pyx_v_dPhir = 0;
  PyObject *__pyx_v_NRPhi = 0;
  PyObject *__pyx_v_margin = 0;
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("_Ves_Vmesh_Phi_numpy (wrapper)", 0);
  {
    static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_DPhi,&__pyx_n_s_dPhir,&__pyx_n_s_NRPhi,&__pyx_n_s_margin,0};
    PyObject* values[4] = {0,0,0,0};
    values[3] = ((PyObject *)__pyx_float_1_eneg_9);
    if (unlikely(__pyx_kwds)) {
      Py_ssize_t kw_args;
      const Py_ssize_t pos_args = PyTuple_GET_SIZE(__pyx_args);
      switch (pos_args) {
        case  4: values[3] = PyTuple_GET_ITEM(__pyx_args, 3);
        case  3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
        case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
        case  1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
        case  0: break;
        default: goto __pyx_L5_argtuple_error;
      }
      kw_args = PyDict_Size(__pyx_kwds);
      switch (pos_args) {
        case  0:
        if (likely((values[0] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_DPhi)) != 0)) kw_args--;
        else goto __pyx_L5_argtuple_error;
        case  1:
        if (likely((values[1] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_dPhir)) != 0)) kw_args--;
        else {
          __Pyx_RaiseArgtupleInvalid("_Ves_Vmesh_Phi_numpy", 0, 3, 4, 1); __PYX_ERR(0, 46, __pyx_L3_error)
        }
        case  2:
        if (likely((values[2] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_NRPhi)) != 0)) kw_args--;
        else {
          __Pyx_RaiseArgtupleInvalid("_Ves_Vmesh_Phi_numpy", 0, 3, 4, 2); __PYX_ERR(0, 46, __pyx_L3_error)
        }
        case  3:
        if (kw_args > 0) {
          PyObject* value = PyDict_GetItem(__pyx_kwds, __pyx_n_s_margin);
          if (value) { values[3] = value; kw_args--; }
        }
      }
      if (unlikely(kw_args > 0)) {
        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "_Ves_Vmesh_Phi_numpy") < 0)) __PYX_ERR(0, 46, __pyx_L3_error)
      }
    } else {
      switch (PyTuple_GET_SIZE(__pyx_args)) {
        case  4: values[3] = PyTuple_GET_ITEM(__pyx_args, 3);
        case  3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);
        values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
        values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
        break;
        default: goto __pyx_L5_argtuple_error;
      }
    }
    __pyx_v_DPhi = values[0];
    __pyx_v_dPhir = values[1];
    __pyx_v_NRPhi = values[2];
    __pyx_v_margin = values[3];
  }
  goto __pyx_L4_argument_unpacking_done;
  __pyx_L5_argtuple_error:;
  __Pyx_RaiseArgtupleInvalid("_Ves_Vmesh_Phi_numpy", 0, 3, 4, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 46, __pyx_L3_error)
  __pyx_L3_error:;
  __Pyx_AddTraceback("_cython_magic_1e9ae0b4a9168851b26e8dfa94f3a771._Ves_Vmesh_Phi_numpy", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __Pyx_RefNannyFinishContext();
  return NULL;
  __pyx_L4_argument_unpacking_done:;
  __pyx_r = __pyx_pf_46_cython_magic_1e9ae0b4a9168851b26e8dfa94f3a771_2_Ves_Vmesh_Phi_numpy(__pyx_self, __pyx_v_DPhi, __pyx_v_dPhir, __pyx_v_NRPhi, __pyx_v_margin);

  /* function exit code */
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_pf_46_cython_magic_1e9ae0b4a9168851b26e8dfa94f3a771_2_Ves_Vmesh_Phi_numpy(CYTHON_UNUSED PyObject *__pyx_self, PyObject *__pyx_v_DPhi, PyObject *__pyx_v_dPhir, PyObject *__pyx_v_NRPhi, PyObject *__pyx_v_margin) {
  PyObject *__pyx_v_nPhi0 = NULL;
  PyObject *__pyx_v_nPhi1 = NULL;
  PyObject *__pyx_v_abs0 = NULL;
  PyObject *__pyx_v_ind0 = NULL;
  PyObject *__pyx_v_abs1 = NULL;
  PyObject *__pyx_v_ind1 = NULL;
  PyObject *__pyx_r = NULL;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("_Ves_Vmesh_Phi_numpy", 0);
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_1);
  __Pyx_XDECREF(__pyx_t_2);
  __Pyx_XDECREF(__pyx_t_3);
  __Pyx_XDECREF(__pyx_t_4);
  __Pyx_XDECREF(__pyx_t_5);
  __Pyx_XDECREF(__pyx_t_6);
  __Pyx_AddTraceback("_cython_magic_1e9ae0b4a9168851b26e8dfa94f3a771._Ves_Vmesh_Phi_numpy", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  __pyx_L0:;
  __Pyx_XDECREF(__pyx_v_nPhi0);
  __Pyx_XDECREF(__pyx_v_nPhi1);
  __Pyx_XDECREF(__pyx_v_abs0);
  __Pyx_XDECREF(__pyx_v_ind0);
  __Pyx_XDECREF(__pyx_v_abs1);
  __Pyx_XDECREF(__pyx_v_ind1);
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
/* … */
  __pyx_tuple__25 = PyTuple_Pack(10, __pyx_n_s_DPhi, __pyx_n_s_dPhir, __pyx_n_s_NRPhi, __pyx_n_s_margin, __pyx_n_s_nPhi0, __pyx_n_s_nPhi1, __pyx_n_s_abs0, __pyx_n_s_ind0, __pyx_n_s_abs1, __pyx_n_s_ind1); if (unlikely(!__pyx_tuple__25)) __PYX_ERR(0, 46, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_tuple__25);
  __Pyx_GIVEREF(__pyx_tuple__25);
/* … */
  __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_46_cython_magic_1e9ae0b4a9168851b26e8dfa94f3a771_3_Ves_Vmesh_Phi_numpy, NULL, __pyx_n_s_cython_magic_1e9ae0b4a9168851b2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 46, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_Ves_Vmesh_Phi_numpy, __pyx_t_1) < 0) __PYX_ERR(0, 46, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_codeobj__26 = (PyObject*)__Pyx_PyCode_New(4, 0, 10, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__25, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_Home_DV226270_cache_ipython_cyt, __pyx_n_s_Ves_Vmesh_Phi_numpy, 46, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__26)) __PYX_ERR(0, 46, __pyx_L1_error)
 47:     """ Get the min and max indices of the relevant Phi for each R (to be finished) """
+48:     nPhi0 = np.nan*np.ones((dPhir.size,))
  __pyx_t_1 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 48, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_nan); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 48, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_3 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 48, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_ones); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 48, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_v_dPhir, __pyx_n_s_size); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 48, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_5 = PyTuple_New(1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 48, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_GIVEREF(__pyx_t_3);
  PyTuple_SET_ITEM(__pyx_t_5, 0, __pyx_t_3);
  __pyx_t_3 = 0;
  __pyx_t_3 = NULL;
  if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_4))) {
    __pyx_t_3 = PyMethod_GET_SELF(__pyx_t_4);
    if (likely(__pyx_t_3)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_4);
      __Pyx_INCREF(__pyx_t_3);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_4, function);
    }
  }
  if (!__pyx_t_3) {
    __pyx_t_1 = __Pyx_PyObject_CallOneArg(__pyx_t_4, __pyx_t_5); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 48, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    __Pyx_GOTREF(__pyx_t_1);
  } else {
    #if CYTHON_FAST_PYCALL
    if (PyFunction_Check(__pyx_t_4)) {
      PyObject *__pyx_temp[2] = {__pyx_t_3, __pyx_t_5};
      __pyx_t_1 = __Pyx_PyFunction_FastCall(__pyx_t_4, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 48, __pyx_L1_error)
      __Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;
      __Pyx_GOTREF(__pyx_t_1);
      __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    } else
    #endif
    #if CYTHON_FAST_PYCCALL
    if (__Pyx_PyFastCFunction_Check(__pyx_t_4)) {
      PyObject *__pyx_temp[2] = {__pyx_t_3, __pyx_t_5};
      __pyx_t_1 = __Pyx_PyCFunction_FastCall(__pyx_t_4, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 48, __pyx_L1_error)
      __Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;
      __Pyx_GOTREF(__pyx_t_1);
      __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    } else
    #endif
    {
      __pyx_t_6 = PyTuple_New(1+1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 48, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_6);
      __Pyx_GIVEREF(__pyx_t_3); PyTuple_SET_ITEM(__pyx_t_6, 0, __pyx_t_3); __pyx_t_3 = NULL;
      __Pyx_GIVEREF(__pyx_t_5);
      PyTuple_SET_ITEM(__pyx_t_6, 0+1, __pyx_t_5);
      __pyx_t_5 = 0;
      __pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_4, __pyx_t_6, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 48, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_1);
      __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
    }
  }
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_t_4 = PyNumber_Multiply(__pyx_t_2, __pyx_t_1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 48, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_v_nPhi0 = __pyx_t_4;
  __pyx_t_4 = 0;
+49:     nPhi1 = np.nan*np.ones((dPhir.size,))
  __pyx_t_4 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 49, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_4, __pyx_n_s_nan); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 49, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_t_2 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 49, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_6 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_ones); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 49, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_6);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_v_dPhir, __pyx_n_s_size); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 49, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_5 = PyTuple_New(1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 49, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_GIVEREF(__pyx_t_2);
  PyTuple_SET_ITEM(__pyx_t_5, 0, __pyx_t_2);
  __pyx_t_2 = 0;
  __pyx_t_2 = NULL;
  if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_6))) {
    __pyx_t_2 = PyMethod_GET_SELF(__pyx_t_6);
    if (likely(__pyx_t_2)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_6);
      __Pyx_INCREF(__pyx_t_2);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_6, function);
    }
  }
  if (!__pyx_t_2) {
    __pyx_t_4 = __Pyx_PyObject_CallOneArg(__pyx_t_6, __pyx_t_5); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 49, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    __Pyx_GOTREF(__pyx_t_4);
  } else {
    #if CYTHON_FAST_PYCALL
    if (PyFunction_Check(__pyx_t_6)) {
      PyObject *__pyx_temp[2] = {__pyx_t_2, __pyx_t_5};
      __pyx_t_4 = __Pyx_PyFunction_FastCall(__pyx_t_6, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 49, __pyx_L1_error)
      __Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
      __Pyx_GOTREF(__pyx_t_4);
      __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    } else
    #endif
    #if CYTHON_FAST_PYCCALL
    if (__Pyx_PyFastCFunction_Check(__pyx_t_6)) {
      PyObject *__pyx_temp[2] = {__pyx_t_2, __pyx_t_5};
      __pyx_t_4 = __Pyx_PyCFunction_FastCall(__pyx_t_6, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 49, __pyx_L1_error)
      __Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
      __Pyx_GOTREF(__pyx_t_4);
      __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    } else
    #endif
    {
      __pyx_t_3 = PyTuple_New(1+1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 49, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_3);
      __Pyx_GIVEREF(__pyx_t_2); PyTuple_SET_ITEM(__pyx_t_3, 0, __pyx_t_2); __pyx_t_2 = NULL;
      __Pyx_GIVEREF(__pyx_t_5);
      PyTuple_SET_ITEM(__pyx_t_3, 0+1, __pyx_t_5);
      __pyx_t_5 = 0;
      __pyx_t_4 = __Pyx_PyObject_Call(__pyx_t_6, __pyx_t_3, NULL); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 49, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_4);
      __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    }
  }
  __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
  __pyx_t_6 = PyNumber_Multiply(__pyx_t_1, __pyx_t_4); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 49, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_6);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_v_nPhi1 = __pyx_t_6;
  __pyx_t_6 = 0;
+50:     abs0 = np.abs(DPhi[0]+np.pi)
  __pyx_t_4 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 50, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_4, __pyx_n_s_abs); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 50, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_t_4 = __Pyx_GetItemInt(__pyx_v_DPhi, 0, long, 1, __Pyx_PyInt_From_long, 0, 0, 1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 50, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_3 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 50, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_pi); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 50, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = PyNumber_Add(__pyx_t_4, __pyx_t_5); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 50, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __pyx_t_5 = NULL;
  if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_1))) {
    __pyx_t_5 = PyMethod_GET_SELF(__pyx_t_1);
    if (likely(__pyx_t_5)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_1);
      __Pyx_INCREF(__pyx_t_5);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_1, function);
    }
  }
  if (!__pyx_t_5) {
    __pyx_t_6 = __Pyx_PyObject_CallOneArg(__pyx_t_1, __pyx_t_3); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 50, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    __Pyx_GOTREF(__pyx_t_6);
  } else {
    #if CYTHON_FAST_PYCALL
    if (PyFunction_Check(__pyx_t_1)) {
      PyObject *__pyx_temp[2] = {__pyx_t_5, __pyx_t_3};
      __pyx_t_6 = __Pyx_PyFunction_FastCall(__pyx_t_1, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 50, __pyx_L1_error)
      __Pyx_XDECREF(__pyx_t_5); __pyx_t_5 = 0;
      __Pyx_GOTREF(__pyx_t_6);
      __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    } else
    #endif
    #if CYTHON_FAST_PYCCALL
    if (__Pyx_PyFastCFunction_Check(__pyx_t_1)) {
      PyObject *__pyx_temp[2] = {__pyx_t_5, __pyx_t_3};
      __pyx_t_6 = __Pyx_PyCFunction_FastCall(__pyx_t_1, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 50, __pyx_L1_error)
      __Pyx_XDECREF(__pyx_t_5); __pyx_t_5 = 0;
      __Pyx_GOTREF(__pyx_t_6);
      __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    } else
    #endif
    {
      __pyx_t_4 = PyTuple_New(1+1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 50, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_4);
      __Pyx_GIVEREF(__pyx_t_5); PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_t_5); __pyx_t_5 = NULL;
      __Pyx_GIVEREF(__pyx_t_3);
      PyTuple_SET_ITEM(__pyx_t_4, 0+1, __pyx_t_3);
      __pyx_t_3 = 0;
      __pyx_t_6 = __Pyx_PyObject_Call(__pyx_t_1, __pyx_t_4, NULL); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 50, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_6);
      __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    }
  }
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_v_abs0 = __pyx_t_6;
  __pyx_t_6 = 0;
+51:     ind0 = abs0-dPhir*np.floor(abs0/dPhir)<margin*dPhir
  __pyx_t_1 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 51, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_floor); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 51, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = __Pyx_PyNumber_Divide(__pyx_v_abs0, __pyx_v_dPhir); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 51, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_3 = NULL;
  if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_4))) {
    __pyx_t_3 = PyMethod_GET_SELF(__pyx_t_4);
    if (likely(__pyx_t_3)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_4);
      __Pyx_INCREF(__pyx_t_3);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_4, function);
    }
  }
  if (!__pyx_t_3) {
    __pyx_t_6 = __Pyx_PyObject_CallOneArg(__pyx_t_4, __pyx_t_1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 51, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
    __Pyx_GOTREF(__pyx_t_6);
  } else {
    #if CYTHON_FAST_PYCALL
    if (PyFunction_Check(__pyx_t_4)) {
      PyObject *__pyx_temp[2] = {__pyx_t_3, __pyx_t_1};
      __pyx_t_6 = __Pyx_PyFunction_FastCall(__pyx_t_4, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 51, __pyx_L1_error)
      __Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;
      __Pyx_GOTREF(__pyx_t_6);
      __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
    } else
    #endif
    #if CYTHON_FAST_PYCCALL
    if (__Pyx_PyFastCFunction_Check(__pyx_t_4)) {
      PyObject *__pyx_temp[2] = {__pyx_t_3, __pyx_t_1};
      __pyx_t_6 = __Pyx_PyCFunction_FastCall(__pyx_t_4, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 51, __pyx_L1_error)
      __Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;
      __Pyx_GOTREF(__pyx_t_6);
      __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
    } else
    #endif
    {
      __pyx_t_5 = PyTuple_New(1+1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 51, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_5);
      __Pyx_GIVEREF(__pyx_t_3); PyTuple_SET_ITEM(__pyx_t_5, 0, __pyx_t_3); __pyx_t_3 = NULL;
      __Pyx_GIVEREF(__pyx_t_1);
      PyTuple_SET_ITEM(__pyx_t_5, 0+1, __pyx_t_1);
      __pyx_t_1 = 0;
      __pyx_t_6 = __Pyx_PyObject_Call(__pyx_t_4, __pyx_t_5, NULL); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 51, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_6);
      __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    }
  }
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_t_4 = PyNumber_Multiply(__pyx_v_dPhir, __pyx_t_6); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 51, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
  __pyx_t_6 = PyNumber_Subtract(__pyx_v_abs0, __pyx_t_4); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 51, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_6);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_t_4 = PyNumber_Multiply(__pyx_v_margin, __pyx_v_dPhir); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 51, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_5 = PyObject_RichCompare(__pyx_t_6, __pyx_t_4, Py_LT); __Pyx_XGOTREF(__pyx_t_5); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 51, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_v_ind0 = __pyx_t_5;
  __pyx_t_5 = 0;
+52:     nPhi0[ind0] = np.round((DPhi[0]+np.pi)/dPhir[ind0])
  __pyx_t_4 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 52, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_6 = __Pyx_PyObject_GetAttrStr(__pyx_t_4, __pyx_n_s_round); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 52, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_6);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_t_4 = __Pyx_GetItemInt(__pyx_v_DPhi, 0, long, 1, __Pyx_PyInt_From_long, 0, 0, 1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 52, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_1 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 52, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_pi); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 52, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = PyNumber_Add(__pyx_t_4, __pyx_t_3); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 52, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = PyObject_GetItem(__pyx_v_dPhir, __pyx_v_ind0); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 52, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_4 = __Pyx_PyNumber_Divide(__pyx_t_1, __pyx_t_3); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 52, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = NULL;
  if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_6))) {
    __pyx_t_3 = PyMethod_GET_SELF(__pyx_t_6);
    if (likely(__pyx_t_3)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_6);
      __Pyx_INCREF(__pyx_t_3);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_6, function);
    }
  }
  if (!__pyx_t_3) {
    __pyx_t_5 = __Pyx_PyObject_CallOneArg(__pyx_t_6, __pyx_t_4); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 52, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    __Pyx_GOTREF(__pyx_t_5);
  } else {
    #if CYTHON_FAST_PYCALL
    if (PyFunction_Check(__pyx_t_6)) {
      PyObject *__pyx_temp[2] = {__pyx_t_3, __pyx_t_4};
      __pyx_t_5 = __Pyx_PyFunction_FastCall(__pyx_t_6, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 52, __pyx_L1_error)
      __Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;
      __Pyx_GOTREF(__pyx_t_5);
      __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    } else
    #endif
    #if CYTHON_FAST_PYCCALL
    if (__Pyx_PyFastCFunction_Check(__pyx_t_6)) {
      PyObject *__pyx_temp[2] = {__pyx_t_3, __pyx_t_4};
      __pyx_t_5 = __Pyx_PyCFunction_FastCall(__pyx_t_6, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 52, __pyx_L1_error)
      __Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;
      __Pyx_GOTREF(__pyx_t_5);
      __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    } else
    #endif
    {
      __pyx_t_1 = PyTuple_New(1+1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 52, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_1);
      __Pyx_GIVEREF(__pyx_t_3); PyTuple_SET_ITEM(__pyx_t_1, 0, __pyx_t_3); __pyx_t_3 = NULL;
      __Pyx_GIVEREF(__pyx_t_4);
      PyTuple_SET_ITEM(__pyx_t_1, 0+1, __pyx_t_4);
      __pyx_t_4 = 0;
      __pyx_t_5 = __Pyx_PyObject_Call(__pyx_t_6, __pyx_t_1, NULL); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 52, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_5);
      __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
    }
  }
  __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
  if (unlikely(PyObject_SetItem(__pyx_v_nPhi0, __pyx_v_ind0, __pyx_t_5) < 0)) __PYX_ERR(0, 52, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
+53:     nPhi0[~ind0] = np.floor((DPhi[0]+np.pi)/dPhir[~ind0])
  __pyx_t_6 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 53, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_6);
  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_6, __pyx_n_s_floor); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 53, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
  __pyx_t_6 = __Pyx_GetItemInt(__pyx_v_DPhi, 0, long, 1, __Pyx_PyInt_From_long, 0, 0, 1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 53, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_6);
  __pyx_t_4 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 53, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_4, __pyx_n_s_pi); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 53, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_t_4 = PyNumber_Add(__pyx_t_6, __pyx_t_3); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 53, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = PyNumber_Invert(__pyx_v_ind0); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 53, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_6 = PyObject_GetItem(__pyx_v_dPhir, __pyx_t_3); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 53, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_6);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = __Pyx_PyNumber_Divide(__pyx_t_4, __pyx_t_6); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 53, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
  __pyx_t_6 = NULL;
  if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_1))) {
    __pyx_t_6 = PyMethod_GET_SELF(__pyx_t_1);
    if (likely(__pyx_t_6)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_1);
      __Pyx_INCREF(__pyx_t_6);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_1, function);
    }
  }
  if (!__pyx_t_6) {
    __pyx_t_5 = __Pyx_PyObject_CallOneArg(__pyx_t_1, __pyx_t_3); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 53, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    __Pyx_GOTREF(__pyx_t_5);
  } else {
    #if CYTHON_FAST_PYCALL
    if (PyFunction_Check(__pyx_t_1)) {
      PyObject *__pyx_temp[2] = {__pyx_t_6, __pyx_t_3};
      __pyx_t_5 = __Pyx_PyFunction_FastCall(__pyx_t_1, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 53, __pyx_L1_error)
      __Pyx_XDECREF(__pyx_t_6); __pyx_t_6 = 0;
      __Pyx_GOTREF(__pyx_t_5);
      __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    } else
    #endif
    #if CYTHON_FAST_PYCCALL
    if (__Pyx_PyFastCFunction_Check(__pyx_t_1)) {
      PyObject *__pyx_temp[2] = {__pyx_t_6, __pyx_t_3};
      __pyx_t_5 = __Pyx_PyCFunction_FastCall(__pyx_t_1, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 53, __pyx_L1_error)
      __Pyx_XDECREF(__pyx_t_6); __pyx_t_6 = 0;
      __Pyx_GOTREF(__pyx_t_5);
      __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    } else
    #endif
    {
      __pyx_t_4 = PyTuple_New(1+1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 53, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_4);
      __Pyx_GIVEREF(__pyx_t_6); PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_t_6); __pyx_t_6 = NULL;
      __Pyx_GIVEREF(__pyx_t_3);
      PyTuple_SET_ITEM(__pyx_t_4, 0+1, __pyx_t_3);
      __pyx_t_3 = 0;
      __pyx_t_5 = __Pyx_PyObject_Call(__pyx_t_1, __pyx_t_4, NULL); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 53, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_5);
      __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    }
  }
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = PyNumber_Invert(__pyx_v_ind0); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 53, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (unlikely(PyObject_SetItem(__pyx_v_nPhi0, __pyx_t_1, __pyx_t_5) < 0)) __PYX_ERR(0, 53, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
+54:     abs1 = np.abs(DPhi[1]+np.pi)
  __pyx_t_1 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 54, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_abs); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 54, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = __Pyx_GetItemInt(__pyx_v_DPhi, 1, long, 1, __Pyx_PyInt_From_long, 0, 0, 1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 54, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_3 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 54, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_6 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_pi); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 54, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_6);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = PyNumber_Add(__pyx_t_1, __pyx_t_6); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 54, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
  __pyx_t_6 = NULL;
  if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_4))) {
    __pyx_t_6 = PyMethod_GET_SELF(__pyx_t_4);
    if (likely(__pyx_t_6)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_4);
      __Pyx_INCREF(__pyx_t_6);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_4, function);
    }
  }
  if (!__pyx_t_6) {
    __pyx_t_5 = __Pyx_PyObject_CallOneArg(__pyx_t_4, __pyx_t_3); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 54, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    __Pyx_GOTREF(__pyx_t_5);
  } else {
    #if CYTHON_FAST_PYCALL
    if (PyFunction_Check(__pyx_t_4)) {
      PyObject *__pyx_temp[2] = {__pyx_t_6, __pyx_t_3};
      __pyx_t_5 = __Pyx_PyFunction_FastCall(__pyx_t_4, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 54, __pyx_L1_error)
      __Pyx_XDECREF(__pyx_t_6); __pyx_t_6 = 0;
      __Pyx_GOTREF(__pyx_t_5);
      __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    } else
    #endif
    #if CYTHON_FAST_PYCCALL
    if (__Pyx_PyFastCFunction_Check(__pyx_t_4)) {
      PyObject *__pyx_temp[2] = {__pyx_t_6, __pyx_t_3};
      __pyx_t_5 = __Pyx_PyCFunction_FastCall(__pyx_t_4, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 54, __pyx_L1_error)
      __Pyx_XDECREF(__pyx_t_6); __pyx_t_6 = 0;
      __Pyx_GOTREF(__pyx_t_5);
      __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    } else
    #endif
    {
      __pyx_t_1 = PyTuple_New(1+1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 54, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_1);
      __Pyx_GIVEREF(__pyx_t_6); PyTuple_SET_ITEM(__pyx_t_1, 0, __pyx_t_6); __pyx_t_6 = NULL;
      __Pyx_GIVEREF(__pyx_t_3);
      PyTuple_SET_ITEM(__pyx_t_1, 0+1, __pyx_t_3);
      __pyx_t_3 = 0;
      __pyx_t_5 = __Pyx_PyObject_Call(__pyx_t_4, __pyx_t_1, NULL); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 54, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_5);
      __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
    }
  }
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_v_abs1 = __pyx_t_5;
  __pyx_t_5 = 0;
+55:     ind1 = abs1-dPhir*np.floor(abs1/dPhir)<margin*dPhir
  __pyx_t_4 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 55, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_4, __pyx_n_s_floor); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 55, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_t_4 = __Pyx_PyNumber_Divide(__pyx_v_abs1, __pyx_v_dPhir); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 55, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_3 = NULL;
  if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_1))) {
    __pyx_t_3 = PyMethod_GET_SELF(__pyx_t_1);
    if (likely(__pyx_t_3)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_1);
      __Pyx_INCREF(__pyx_t_3);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_1, function);
    }
  }
  if (!__pyx_t_3) {
    __pyx_t_5 = __Pyx_PyObject_CallOneArg(__pyx_t_1, __pyx_t_4); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 55, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    __Pyx_GOTREF(__pyx_t_5);
  } else {
    #if CYTHON_FAST_PYCALL
    if (PyFunction_Check(__pyx_t_1)) {
      PyObject *__pyx_temp[2] = {__pyx_t_3, __pyx_t_4};
      __pyx_t_5 = __Pyx_PyFunction_FastCall(__pyx_t_1, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 55, __pyx_L1_error)
      __Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;
      __Pyx_GOTREF(__pyx_t_5);
      __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    } else
    #endif
    #if CYTHON_FAST_PYCCALL
    if (__Pyx_PyFastCFunction_Check(__pyx_t_1)) {
      PyObject *__pyx_temp[2] = {__pyx_t_3, __pyx_t_4};
      __pyx_t_5 = __Pyx_PyCFunction_FastCall(__pyx_t_1, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 55, __pyx_L1_error)
      __Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;
      __Pyx_GOTREF(__pyx_t_5);
      __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    } else
    #endif
    {
      __pyx_t_6 = PyTuple_New(1+1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 55, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_6);
      __Pyx_GIVEREF(__pyx_t_3); PyTuple_SET_ITEM(__pyx_t_6, 0, __pyx_t_3); __pyx_t_3 = NULL;
      __Pyx_GIVEREF(__pyx_t_4);
      PyTuple_SET_ITEM(__pyx_t_6, 0+1, __pyx_t_4);
      __pyx_t_4 = 0;
      __pyx_t_5 = __Pyx_PyObject_Call(__pyx_t_1, __pyx_t_6, NULL); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 55, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_5);
      __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
    }
  }
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = PyNumber_Multiply(__pyx_v_dPhir, __pyx_t_5); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 55, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __pyx_t_5 = PyNumber_Subtract(__pyx_v_abs1, __pyx_t_1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 55, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = PyNumber_Multiply(__pyx_v_margin, __pyx_v_dPhir); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 55, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_6 = PyObject_RichCompare(__pyx_t_5, __pyx_t_1, Py_LT); __Pyx_XGOTREF(__pyx_t_6); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 55, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_v_ind1 = __pyx_t_6;
  __pyx_t_6 = 0;
+56:     nPhi1[ind1] = np.round((DPhi[1]+np.pi)/dPhir[ind1])-1
  __pyx_t_1 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 56, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_round); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 56, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = __Pyx_GetItemInt(__pyx_v_DPhi, 1, long, 1, __Pyx_PyInt_From_long, 0, 0, 1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 56, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_4 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 56, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_4, __pyx_n_s_pi); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 56, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_t_4 = PyNumber_Add(__pyx_t_1, __pyx_t_3); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 56, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = PyObject_GetItem(__pyx_v_dPhir, __pyx_v_ind1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 56, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_1 = __Pyx_PyNumber_Divide(__pyx_t_4, __pyx_t_3); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 56, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = NULL;
  if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_5))) {
    __pyx_t_3 = PyMethod_GET_SELF(__pyx_t_5);
    if (likely(__pyx_t_3)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_5);
      __Pyx_INCREF(__pyx_t_3);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_5, function);
    }
  }
  if (!__pyx_t_3) {
    __pyx_t_6 = __Pyx_PyObject_CallOneArg(__pyx_t_5, __pyx_t_1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 56, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
    __Pyx_GOTREF(__pyx_t_6);
  } else {
    #if CYTHON_FAST_PYCALL
    if (PyFunction_Check(__pyx_t_5)) {
      PyObject *__pyx_temp[2] = {__pyx_t_3, __pyx_t_1};
      __pyx_t_6 = __Pyx_PyFunction_FastCall(__pyx_t_5, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 56, __pyx_L1_error)
      __Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;
      __Pyx_GOTREF(__pyx_t_6);
      __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
    } else
    #endif
    #if CYTHON_FAST_PYCCALL
    if (__Pyx_PyFastCFunction_Check(__pyx_t_5)) {
      PyObject *__pyx_temp[2] = {__pyx_t_3, __pyx_t_1};
      __pyx_t_6 = __Pyx_PyCFunction_FastCall(__pyx_t_5, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 56, __pyx_L1_error)
      __Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;
      __Pyx_GOTREF(__pyx_t_6);
      __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
    } else
    #endif
    {
      __pyx_t_4 = PyTuple_New(1+1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 56, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_4);
      __Pyx_GIVEREF(__pyx_t_3); PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_t_3); __pyx_t_3 = NULL;
      __Pyx_GIVEREF(__pyx_t_1);
      PyTuple_SET_ITEM(__pyx_t_4, 0+1, __pyx_t_1);
      __pyx_t_1 = 0;
      __pyx_t_6 = __Pyx_PyObject_Call(__pyx_t_5, __pyx_t_4, NULL); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 56, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_6);
      __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    }
  }
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __pyx_t_5 = __Pyx_PyInt_SubtractObjC(__pyx_t_6, __pyx_int_1, 1, 0); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 56, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
  if (unlikely(PyObject_SetItem(__pyx_v_nPhi1, __pyx_v_ind1, __pyx_t_5) < 0)) __PYX_ERR(0, 56, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
+57:     nPhi1[~ind1] = np.floor((DPhi[1]+np.pi)/dPhir[~ind1])
  __pyx_t_6 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 57, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_6);
  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_6, __pyx_n_s_floor); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 57, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
  __pyx_t_6 = __Pyx_GetItemInt(__pyx_v_DPhi, 1, long, 1, __Pyx_PyInt_From_long, 0, 0, 1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 57, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_6);
  __pyx_t_1 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 57, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_pi); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 57, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = PyNumber_Add(__pyx_t_6, __pyx_t_3); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 57, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = PyNumber_Invert(__pyx_v_ind1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 57, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_6 = PyObject_GetItem(__pyx_v_dPhir, __pyx_t_3); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 57, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_6);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = __Pyx_PyNumber_Divide(__pyx_t_1, __pyx_t_6); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 57, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
  __pyx_t_6 = NULL;
  if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_4))) {
    __pyx_t_6 = PyMethod_GET_SELF(__pyx_t_4);
    if (likely(__pyx_t_6)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_4);
      __Pyx_INCREF(__pyx_t_6);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_4, function);
    }
  }
  if (!__pyx_t_6) {
    __pyx_t_5 = __Pyx_PyObject_CallOneArg(__pyx_t_4, __pyx_t_3); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 57, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    __Pyx_GOTREF(__pyx_t_5);
  } else {
    #if CYTHON_FAST_PYCALL
    if (PyFunction_Check(__pyx_t_4)) {
      PyObject *__pyx_temp[2] = {__pyx_t_6, __pyx_t_3};
      __pyx_t_5 = __Pyx_PyFunction_FastCall(__pyx_t_4, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 57, __pyx_L1_error)
      __Pyx_XDECREF(__pyx_t_6); __pyx_t_6 = 0;
      __Pyx_GOTREF(__pyx_t_5);
      __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    } else
    #endif
    #if CYTHON_FAST_PYCCALL
    if (__Pyx_PyFastCFunction_Check(__pyx_t_4)) {
      PyObject *__pyx_temp[2] = {__pyx_t_6, __pyx_t_3};
      __pyx_t_5 = __Pyx_PyCFunction_FastCall(__pyx_t_4, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 57, __pyx_L1_error)
      __Pyx_XDECREF(__pyx_t_6); __pyx_t_6 = 0;
      __Pyx_GOTREF(__pyx_t_5);
      __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    } else
    #endif
    {
      __pyx_t_1 = PyTuple_New(1+1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 57, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_1);
      __Pyx_GIVEREF(__pyx_t_6); PyTuple_SET_ITEM(__pyx_t_1, 0, __pyx_t_6); __pyx_t_6 = NULL;
      __Pyx_GIVEREF(__pyx_t_3);
      PyTuple_SET_ITEM(__pyx_t_1, 0+1, __pyx_t_3);
      __pyx_t_3 = 0;
      __pyx_t_5 = __Pyx_PyObject_Call(__pyx_t_4, __pyx_t_1, NULL); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 57, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_5);
      __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
    }
  }
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_t_4 = PyNumber_Invert(__pyx_v_ind1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 57, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  if (unlikely(PyObject_SetItem(__pyx_v_nPhi1, __pyx_t_4, __pyx_t_5) < 0)) __PYX_ERR(0, 57, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
+58:     assert np.all(nPhi0>=0) and np.all(nPhi0<=NRPhi)
  #ifndef CYTHON_WITHOUT_ASSERTIONS
  if (unlikely(!Py_OptimizeFlag)) {
    __pyx_t_4 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 58, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_4);
    __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_4, __pyx_n_s_all); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 58, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    __pyx_t_4 = PyObject_RichCompare(__pyx_v_nPhi0, __pyx_int_0, Py_GE); __Pyx_XGOTREF(__pyx_t_4); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 58, __pyx_L1_error)
    __pyx_t_3 = NULL;
    if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_1))) {
      __pyx_t_3 = PyMethod_GET_SELF(__pyx_t_1);
      if (likely(__pyx_t_3)) {
        PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_1);
        __Pyx_INCREF(__pyx_t_3);
        __Pyx_INCREF(function);
        __Pyx_DECREF_SET(__pyx_t_1, function);
      }
    }
    if (!__pyx_t_3) {
      __pyx_t_5 = __Pyx_PyObject_CallOneArg(__pyx_t_1, __pyx_t_4); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 58, __pyx_L1_error)
      __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
      __Pyx_GOTREF(__pyx_t_5);
    } else {
      #if CYTHON_FAST_PYCALL
      if (PyFunction_Check(__pyx_t_1)) {
        PyObject *__pyx_temp[2] = {__pyx_t_3, __pyx_t_4};
        __pyx_t_5 = __Pyx_PyFunction_FastCall(__pyx_t_1, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 58, __pyx_L1_error)
        __Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;
        __Pyx_GOTREF(__pyx_t_5);
        __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
      } else
      #endif
      #if CYTHON_FAST_PYCCALL
      if (__Pyx_PyFastCFunction_Check(__pyx_t_1)) {
        PyObject *__pyx_temp[2] = {__pyx_t_3, __pyx_t_4};
        __pyx_t_5 = __Pyx_PyCFunction_FastCall(__pyx_t_1, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 58, __pyx_L1_error)
        __Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;
        __Pyx_GOTREF(__pyx_t_5);
        __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
      } else
      #endif
      {
        __pyx_t_6 = PyTuple_New(1+1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 58, __pyx_L1_error)
        __Pyx_GOTREF(__pyx_t_6);
        __Pyx_GIVEREF(__pyx_t_3); PyTuple_SET_ITEM(__pyx_t_6, 0, __pyx_t_3); __pyx_t_3 = NULL;
        __Pyx_GIVEREF(__pyx_t_4);
        PyTuple_SET_ITEM(__pyx_t_6, 0+1, __pyx_t_4);
        __pyx_t_4 = 0;
        __pyx_t_5 = __Pyx_PyObject_Call(__pyx_t_1, __pyx_t_6, NULL); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 58, __pyx_L1_error)
        __Pyx_GOTREF(__pyx_t_5);
        __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
      }
    }
    __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
    __pyx_t_8 = __Pyx_PyObject_IsTrue(__pyx_t_5); if (unlikely(__pyx_t_8 < 0)) __PYX_ERR(0, 58, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    if (__pyx_t_8) {
    } else {
      __pyx_t_7 = __pyx_t_8;
      goto __pyx_L3_bool_binop_done;
    }
    __pyx_t_1 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 58, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    __pyx_t_6 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_all); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 58, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_6);
    __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
    __pyx_t_1 = PyObject_RichCompare(__pyx_v_nPhi0, __pyx_v_NRPhi, Py_LE); __Pyx_XGOTREF(__pyx_t_1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 58, __pyx_L1_error)
    __pyx_t_4 = NULL;
    if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_6))) {
      __pyx_t_4 = PyMethod_GET_SELF(__pyx_t_6);
      if (likely(__pyx_t_4)) {
        PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_6);
        __Pyx_INCREF(__pyx_t_4);
        __Pyx_INCREF(function);
        __Pyx_DECREF_SET(__pyx_t_6, function);
      }
    }
    if (!__pyx_t_4) {
      __pyx_t_5 = __Pyx_PyObject_CallOneArg(__pyx_t_6, __pyx_t_1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 58, __pyx_L1_error)
      __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
      __Pyx_GOTREF(__pyx_t_5);
    } else {
      #if CYTHON_FAST_PYCALL
      if (PyFunction_Check(__pyx_t_6)) {
        PyObject *__pyx_temp[2] = {__pyx_t_4, __pyx_t_1};
        __pyx_t_5 = __Pyx_PyFunction_FastCall(__pyx_t_6, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 58, __pyx_L1_error)
        __Pyx_XDECREF(__pyx_t_4); __pyx_t_4 = 0;
        __Pyx_GOTREF(__pyx_t_5);
        __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
      } else
      #endif
      #if CYTHON_FAST_PYCCALL
      if (__Pyx_PyFastCFunction_Check(__pyx_t_6)) {
        PyObject *__pyx_temp[2] = {__pyx_t_4, __pyx_t_1};
        __pyx_t_5 = __Pyx_PyCFunction_FastCall(__pyx_t_6, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 58, __pyx_L1_error)
        __Pyx_XDECREF(__pyx_t_4); __pyx_t_4 = 0;
        __Pyx_GOTREF(__pyx_t_5);
        __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
      } else
      #endif
      {
        __pyx_t_3 = PyTuple_New(1+1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 58, __pyx_L1_error)
        __Pyx_GOTREF(__pyx_t_3);
        __Pyx_GIVEREF(__pyx_t_4); PyTuple_SET_ITEM(__pyx_t_3, 0, __pyx_t_4); __pyx_t_4 = NULL;
        __Pyx_GIVEREF(__pyx_t_1);
        PyTuple_SET_ITEM(__pyx_t_3, 0+1, __pyx_t_1);
        __pyx_t_1 = 0;
        __pyx_t_5 = __Pyx_PyObject_Call(__pyx_t_6, __pyx_t_3, NULL); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 58, __pyx_L1_error)
        __Pyx_GOTREF(__pyx_t_5);
        __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
      }
    }
    __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
    __pyx_t_8 = __Pyx_PyObject_IsTrue(__pyx_t_5); if (unlikely(__pyx_t_8 < 0)) __PYX_ERR(0, 58, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    __pyx_t_7 = __pyx_t_8;
    __pyx_L3_bool_binop_done:;
    if (unlikely(!__pyx_t_7)) {
      PyErr_SetNone(PyExc_AssertionError);
      __PYX_ERR(0, 58, __pyx_L1_error)
    }
  }
  #endif
+59:     assert np.all(nPhi1>=0) and np.all(nPhi1<=NRPhi)
  #ifndef CYTHON_WITHOUT_ASSERTIONS
  if (unlikely(!Py_OptimizeFlag)) {
    __pyx_t_6 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 59, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_6);
    __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_6, __pyx_n_s_all); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 59, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_3);
    __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
    __pyx_t_6 = PyObject_RichCompare(__pyx_v_nPhi1, __pyx_int_0, Py_GE); __Pyx_XGOTREF(__pyx_t_6); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 59, __pyx_L1_error)
    __pyx_t_1 = NULL;
    if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_3))) {
      __pyx_t_1 = PyMethod_GET_SELF(__pyx_t_3);
      if (likely(__pyx_t_1)) {
        PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_3);
        __Pyx_INCREF(__pyx_t_1);
        __Pyx_INCREF(function);
        __Pyx_DECREF_SET(__pyx_t_3, function);
      }
    }
    if (!__pyx_t_1) {
      __pyx_t_5 = __Pyx_PyObject_CallOneArg(__pyx_t_3, __pyx_t_6); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 59, __pyx_L1_error)
      __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
      __Pyx_GOTREF(__pyx_t_5);
    } else {
      #if CYTHON_FAST_PYCALL
      if (PyFunction_Check(__pyx_t_3)) {
        PyObject *__pyx_temp[2] = {__pyx_t_1, __pyx_t_6};
        __pyx_t_5 = __Pyx_PyFunction_FastCall(__pyx_t_3, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 59, __pyx_L1_error)
        __Pyx_XDECREF(__pyx_t_1); __pyx_t_1 = 0;
        __Pyx_GOTREF(__pyx_t_5);
        __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
      } else
      #endif
      #if CYTHON_FAST_PYCCALL
      if (__Pyx_PyFastCFunction_Check(__pyx_t_3)) {
        PyObject *__pyx_temp[2] = {__pyx_t_1, __pyx_t_6};
        __pyx_t_5 = __Pyx_PyCFunction_FastCall(__pyx_t_3, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 59, __pyx_L1_error)
        __Pyx_XDECREF(__pyx_t_1); __pyx_t_1 = 0;
        __Pyx_GOTREF(__pyx_t_5);
        __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
      } else
      #endif
      {
        __pyx_t_4 = PyTuple_New(1+1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 59, __pyx_L1_error)
        __Pyx_GOTREF(__pyx_t_4);
        __Pyx_GIVEREF(__pyx_t_1); PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_t_1); __pyx_t_1 = NULL;
        __Pyx_GIVEREF(__pyx_t_6);
        PyTuple_SET_ITEM(__pyx_t_4, 0+1, __pyx_t_6);
        __pyx_t_6 = 0;
        __pyx_t_5 = __Pyx_PyObject_Call(__pyx_t_3, __pyx_t_4, NULL); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 59, __pyx_L1_error)
        __Pyx_GOTREF(__pyx_t_5);
        __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
      }
    }
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    __pyx_t_8 = __Pyx_PyObject_IsTrue(__pyx_t_5); if (unlikely(__pyx_t_8 < 0)) __PYX_ERR(0, 59, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    if (__pyx_t_8) {
    } else {
      __pyx_t_7 = __pyx_t_8;
      goto __pyx_L5_bool_binop_done;
    }
    __pyx_t_3 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 59, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_3);
    __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_all); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 59, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_4);
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    __pyx_t_3 = PyObject_RichCompare(__pyx_v_nPhi1, __pyx_v_NRPhi, Py_LE); __Pyx_XGOTREF(__pyx_t_3); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 59, __pyx_L1_error)
    __pyx_t_6 = NULL;
    if (CYTHON_UNPACK_METHODS && unlikely(PyMethod_Check(__pyx_t_4))) {
      __pyx_t_6 = PyMethod_GET_SELF(__pyx_t_4);
      if (likely(__pyx_t_6)) {
        PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_4);
        __Pyx_INCREF(__pyx_t_6);
        __Pyx_INCREF(function);
        __Pyx_DECREF_SET(__pyx_t_4, function);
      }
    }
    if (!__pyx_t_6) {
      __pyx_t_5 = __Pyx_PyObject_CallOneArg(__pyx_t_4, __pyx_t_3); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 59, __pyx_L1_error)
      __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
      __Pyx_GOTREF(__pyx_t_5);
    } else {
      #if CYTHON_FAST_PYCALL
      if (PyFunction_Check(__pyx_t_4)) {
        PyObject *__pyx_temp[2] = {__pyx_t_6, __pyx_t_3};
        __pyx_t_5 = __Pyx_PyFunction_FastCall(__pyx_t_4, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 59, __pyx_L1_error)
        __Pyx_XDECREF(__pyx_t_6); __pyx_t_6 = 0;
        __Pyx_GOTREF(__pyx_t_5);
        __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
      } else
      #endif
      #if CYTHON_FAST_PYCCALL
      if (__Pyx_PyFastCFunction_Check(__pyx_t_4)) {
        PyObject *__pyx_temp[2] = {__pyx_t_6, __pyx_t_3};
        __pyx_t_5 = __Pyx_PyCFunction_FastCall(__pyx_t_4, __pyx_temp+1-1, 1+1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 59, __pyx_L1_error)
        __Pyx_XDECREF(__pyx_t_6); __pyx_t_6 = 0;
        __Pyx_GOTREF(__pyx_t_5);
        __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
      } else
      #endif
      {
        __pyx_t_1 = PyTuple_New(1+1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 59, __pyx_L1_error)
        __Pyx_GOTREF(__pyx_t_1);
        __Pyx_GIVEREF(__pyx_t_6); PyTuple_SET_ITEM(__pyx_t_1, 0, __pyx_t_6); __pyx_t_6 = NULL;
        __Pyx_GIVEREF(__pyx_t_3);
        PyTuple_SET_ITEM(__pyx_t_1, 0+1, __pyx_t_3);
        __pyx_t_3 = 0;
        __pyx_t_5 = __Pyx_PyObject_Call(__pyx_t_4, __pyx_t_1, NULL); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 59, __pyx_L1_error)
        __Pyx_GOTREF(__pyx_t_5);
        __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
      }
    }
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    __pyx_t_8 = __Pyx_PyObject_IsTrue(__pyx_t_5); if (unlikely(__pyx_t_8 < 0)) __PYX_ERR(0, 59, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    __pyx_t_7 = __pyx_t_8;
    __pyx_L5_bool_binop_done:;
    if (unlikely(!__pyx_t_7)) {
      PyErr_SetNone(PyExc_AssertionError);
      __PYX_ERR(0, 59, __pyx_L1_error)
    }
  }
  #endif
+60:     return nPhi0, nPhi1
  __Pyx_XDECREF(__pyx_r);
  __pyx_t_5 = PyTuple_New(2); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 60, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_INCREF(__pyx_v_nPhi0);
  __Pyx_GIVEREF(__pyx_v_nPhi0);
  PyTuple_SET_ITEM(__pyx_t_5, 0, __pyx_v_nPhi0);
  __Pyx_INCREF(__pyx_v_nPhi1);
  __Pyx_GIVEREF(__pyx_v_nPhi1);
  PyTuple_SET_ITEM(__pyx_t_5, 1, __pyx_v_nPhi1);
  __pyx_r = __pyx_t_5;
  __pyx_t_5 = 0;
  goto __pyx_L0;



In [ ]: